/*
  Copyright (C) 2011, Kyungjoo Kim
  
  This file is part of UHM (Unassembled HyperMatrix) sparse direct solver.
  
  All rights reserved.

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.
  
  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.
  
  You should have received a copy of the GNU Lesser General Public
  License along with this library; if not, write to the Free Software
  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  Also add information on how to contact you by electronic and paper mail.
  
  
  Kyungjoo Kim
  iamkyungjoo@gmail.com
*/
#include "uhm/common.hxx"
#include "uhm/const.hxx"

#include "uhm/matrix/uhm/matrix.hxx"
#include "uhm/matrix/uhm/fla.hxx"

namespace uhm {
  static int ldl_piv_flat   ( int fs, int ss,
                              linal::Flat_ ATL, linal::Flat_ ATR,
                              linal::Flat_ ABL, linal::Flat_ ABR,
                              linal::Flat_ p );
  
  static int ldl_piv_hier   ( int fs, int ss, 
                              linal::Hier_ ATL, linal::Hier_ ATR,
                              linal::Hier_ ABL, linal::Hier_ ABR,
                              linal::Hier_ p,
                              // ********************************
                              linal::Flat_ FTL, linal::Flat_ FTR,
                              linal::Flat_ FBL );
  
  static int ldl_incpiv_hier( int fs, int ss, 
                              linal::Hier_ ATL, linal::Hier_ ATR,
                              linal::Hier_ ABL, linal::Hier_ ABR,
                              linal::Hier_ p,
                              // ********************************
                              linal::Flat_ FTL, linal::Flat_ FTR,
                              linal::Flat_ FBL );
  
  void Matrix_FLA_::ldl_piv() {

#ifdef UHM_HIER_MATRIX_ENABLE
    // ----------------------------------------------------------
    // ** Hier-Matrix
    // ----------------------------------------------------------
    ldl_piv_hier( this->fs, this->ss, 
                  this->hier.ATL, this->hier.ATR,
                  this->hier.ABL, this->hier.ABR,
                  this->hier.p,
                  this->flat.ATL, this->flat.ATR,
                  this->flat.ABL );
#else
    // ----------------------------------------------------------
    // ** Flat-Matrix - Level Matrix
    // ----------------------------------------------------------
    ldl_piv_flat( this->fs, this->ss, 
                  this->flat.ATL, this->flat.ATR,
                  this->flat.ABL, this->flat.ABR,
                  this->flat.p );
#endif
  }
  
  void Matrix_FLA_::ldl_incpiv() {

#ifdef UHM_HIER_MATRIX_ENABLE
    // ----------------------------------------------------------
    // ** Hier-Matrix
    // ----------------------------------------------------------
    ldl_incpiv_hier( this->fs, this->ss, 
                     this->hier.ATL, this->hier.ATR,
                     this->hier.ABL, this->hier.ABR,
                     this->hier.p,
                     this->flat.ATL, this->flat.ATR,
                     this->flat.ABL );
#else
    fprintf(stderr, "incremental pivoting is not supported in flat matrix\n");
    abort();
#endif
  }
  
  static inline int ldl_piv_flat( int fs, int ss, 
                                  linal::Flat_ ATL, linal::Flat_ ATR,
                                  linal::Flat_ ABL, linal::Flat_ ABR,
                                  linal::Flat_ p ) {
    // ** BLAS 2 : do not expect high performance
    //    If possible, use incremental pivoting strategy   
    if (fs) 
      linal::host::internal::ldl_piv( FLA_LOWER_TRIANGULAR, ~ATL, ~p );
  
    if (fs && ss) {

#pragma omp task firstprivate( ATL )
      {
        FLA_Hermitianize( FLA_LOWER_TRIANGULAR, ~ATL );
        linal::host::internal::apply_diag_scal( FLA_LEFT, FLA_NO_CONJUGATE,
                                                FLA_UPPER_TRIANGULAR, ~ATL );
      }

#pragma omp task firstprivate( p, ABL )
      linal::host::internal::apply_pivots( FLA_RIGHT, FLA_TRANSPOSE, ~p, ~ABL );

#pragma omp taskwait

      FLA_Trsm( FLA_RIGHT, FLA_UPPER_TRIANGULAR,
                FLA_NO_TRANSPOSE, FLA_NONUNIT_DIAG,
                FLA_ONE, ~ATL, ~ABL );

      ATR.create( ABL.get_datatype(), ABL.get_n(), ABL.get_m() );

      FLA_Copyt( FLA_TRANSPOSE, ~ABL, ~ATR );
      linal::host::internal::apply_diag_scal( FLA_LEFT, FLA_NO_CONJUGATE, 
                                              ~ATL, ~ATR );

#pragma omp taskwait

      FLA_Gemm( FLA_NO_TRANSPOSE, FLA_NO_TRANSPOSE,
                FLA_MINUS_ONE, ~ABL, ~ATR, FLA_ONE, ~ABR );
      
      ATR.free();
    }

    return true;
  }


  static inline int ldl_piv_hier( int fs, int ss, 
                                  linal::Hier_ ATL, linal::Hier_ ATR,
                                  linal::Hier_ ABL, linal::Hier_ ABR,
                                  linal::Hier_ p,
                                  linal::Flat_ FTL, linal::Flat_ FTR,
                                  linal::Flat_ FBL ) {
    // ** BLAS 2 : do not expect high performance
    //    If possible, use incremental pivoting strategy   
    if (fs) 
      linal::host::internal::ldl_piv( FLA_LOWER_TRIANGULAR, ~FTL, ~(p.flat()) );
  
    if (fs && ss) {

#pragma omp task firstprivate( FTL )
      {
        FLA_Hermitianize( FLA_LOWER_TRIANGULAR, ~FTL );
        linal::host::internal::apply_diag_scal( FLA_LEFT, FLA_NO_CONJUGATE, 
                                                FLA_UPPER_TRIANGULAR, ~FTL );
      }

#pragma omp task firstprivate( p, ABL )
      linal::apply_pivots( FLA_RIGHT, FLA_TRANSPOSE, p, ABL );

#pragma omp taskwait

      linal::trsm( FLA_RIGHT, FLA_UPPER_TRIANGULAR,
                          FLA_NO_TRANSPOSE, FLA_NONUNIT_DIAG,
                          FLA_ONE, ATL, ABL );

      int bmn = get_hier_blocksize();

      FTR.create( FBL.get_datatype(), FBL.get_n(), FBL.get_m() );
      ATR.create( FTR, bmn, bmn );

      FLA_Copyt( FLA_TRANSPOSE, ~FBL, ~FTR );
      linal::host::internal::apply_diag_scal( FLA_LEFT, FLA_NO_CONJUGATE, 
                                              ~FTL, ~FTR );

#pragma omp taskwait

      linal::gemm_l( FLA_NO_TRANSPOSE, FLA_NO_TRANSPOSE,
                            FLA_MINUS_ONE, ABL, ATR, FLA_ONE, ABR );
      
      ATR.free();
      FTR.free();
    }

    return true;
  }

  static inline int ldl_incpiv_hier( int fs, int ss, 
                                     linal::Hier_ ATL, linal::Hier_ ATR,
                                     linal::Hier_ ABL, linal::Hier_ ABR,
                                     linal::Hier_ p,
                                     linal::Flat_ FTL, linal::Flat_ FTR,
                                     linal::Flat_ FBL ) {
    if (fs) 
      linal::ldl_incpiv( FLA_LOWER_TRIANGULAR, ATL, p );
  
    if (fs && ss) {

#pragma omp task firstprivate( FTL )
      {
        FLA_Hermitianize( FLA_LOWER_TRIANGULAR, ~FTL );
        linal::host::internal::apply_diag_scal( FLA_LEFT, FLA_NO_CONJUGATE,
                                                FLA_UPPER_TRIANGULAR, ~FTL );
      }

#pragma omp task firstprivate( p, ABL )
      linal::apply_pivots( FLA_RIGHT, FLA_TRANSPOSE, p, ABL );

#pragma omp taskwait

      linal::trsm( FLA_RIGHT, FLA_UPPER_TRIANGULAR,
                          FLA_NO_TRANSPOSE, FLA_NONUNIT_DIAG,
                          FLA_ONE, ATL, ABL );

      int bmn = get_hier_blocksize();

      FTR.create( FBL.get_datatype(), FBL.get_n(), FBL.get_m() );
      ATR.create( FTR, bmn, bmn );

      FLA_Copyt( FLA_TRANSPOSE, ~FBL, ~FTR );
      linal::host::internal::apply_diag_scal( FLA_LEFT, FLA_NO_CONJUGATE,
                                              ~FTL, ~FTR );

#pragma omp taskwait

      linal::gemm_l( FLA_NO_TRANSPOSE, FLA_NO_TRANSPOSE,
                            FLA_MINUS_ONE, ABL, ATR, FLA_ONE, ABR );
      
      ATR.free();
      FTR.free();
    }

    return true;
  }

}
